AV  650!?^- 


MEMORANDUM 

RM-5244-PR 

mCU  1987 


MATHEMATICAL  ANALYSIS 
AND  DIGITAL  SIMULATION  OF  THE 
RESP1 RATORY  CONTROL  SYSTEM 

Fred  S.  Goodins,  M.D.,  Jur  e  Buell  and  Alex  J.  Bart 


This  research  is  supported  by  the  United  States  Air  Force  under  Project  R.  ND— Con¬ 
tract  No.  F44620-67-C-0045— monitored  by  the  Directorate  cf  Operational  Requirements 
and  Development  Plans.  Deputy  Chief  of  StafT.  Rt  earch  and  Development,  Hq  USAF. 
Views  or  conclusions  contained  in  this  Memorandum  should  not  be  interpreted  as 
representing  the  official  opinion  or  policy  of  she  United  States  Air  Foret . 


DISTRIBUTION  STATEMENT 
Distribution  of  this  document  is  unlimited. 


WiflD 


(fn/trtdffat 


tANtA  MONICA  •  CAltPGANI  A  •  Hill. 


I  PAS  M 


M 


-iii- 

PREPACE 

This  Memorandum  describes  a  mathematical  model  of  the  human 
respiratory  control  system.  It  was  supported  in  part  by  the  National 
Institutes  of  Health  under  the  following  giants:  GM-09608-04; 

4-K6-HE  14,  187;  and  5-ROI-KE  01626. 

Two  of  the  authors  are  now  at  the  Northwestern  University 
Medical  School:  F.  S.  Grodins,  M.D.,  Ph.D.,  consultant  to  The  RAND 
Corporation,  is  a  Professor  of  Physiology,  and  A.  J.  Bart  is  an 
American  Cancer  Society  Fellow  there. 

A  version  of  this  Memorandum  was  accepted  for  publication  ir 
the  Journal  of  Applied  Physiology. 
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This  Memorandum  expresses  the  basic  material  balance  relationships 
for  the  lung-blood-tissue  gas  transport  and  exchange  system  in  a  set 
of  differential-difference  equations  containing  a  number  of  dependent 
time  delays.  Additional  equations  define  the  chemical  details  of 
transport  and  acid— base  buffering,  concentration  equilibria,  and  blood 


flow  behavior.  Finally,  a  control  function  is  included  defining  the 

dependence  of  ventilation  upon  CSF  (flf*),  and  arterial  (H*)  and  Pn 

u2 

at  the  carotid  chemoceptors .  A  Foreran  program  was  written  for 


convenient  digital  simulation  of  the  responses  of  the  sys.tem  to  a 


wide  variety  of  forcings,  including  CC^  inhalation,  hypoxia  at  sea 
level,  altitude  hypoxia,  and  metabolic  disturbances  in  acid-base 


balance.  Both  dynamic  and  steady-state  behavior  of  the  model  were 


reasonably  realistic. 
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MATHEMATICAL  ANALYSIS  AND  DIGITAL  SIMULATION 
OF  THE  RESPIRATORY  CONTROL  SYSTEM 


!_. _ INTRODUCTION 

Although  recognition  of  the  basic  "closed-loop"  nature  of  the 
respiratory  control  system  was  certainly  implicit  in  the  writings  of 
Haldane  and  Pvt  or*-i »cy  [1],  the  first  explicit  quantitative  formula¬ 
tion  was  made  by  Gray  in  1945  [2]  His  algebraic  model  was  restricted 
to  the  steady-state  responses  of  the  system  to  COp  inhalation, 
arterial  anoxemia,  and  metabolic  disturbances  in  acid-base  balance. 

The  first  dynamic  analysis  of  this  system  appeared  in  1954  [3],  and 
although  it  represented  a  step  forward  in  considering  dynamics,  it 
also  represented  a  step  backward  in  that  it  could  accept  only  a 
single  forcing,  i.e.,  C0£  inhalation.  It  also  represented  an  over¬ 
simplified  treatment  that  nevertheless  led  to  nonlinear  differential 
equations.  At  the  time  of  its  formulation,  the  state  of  the  computer 
art  was  relatively  primitive,  and  exploration  of  a  more  general  and 
realistic  model  was  not  practical. 

With  the  development  of  large  computing  facilities,  interest 
in  the  realistic  dynamic  analysis  of  complex  multivariate  biological 
control  systems  has  increased  rapidly  [4].  Beginning  with  Defares, 
Derkson,  and  Duyff  in  1960  [5],  the  original  dynamic  model  of  the 
respiratory  system  has  been  refined  and  extended  by  several  workers 
using  both  analog  and  digital  simulation  [6-9 3 •  The  present  analysis 
represents  a  further  step  in  this  continuing  study.  The  model  has 
been  made  sufficiently  general  to  acconaaodate  a  variety  of  forcings 
( C(>2  inhalation,  hypoxia  at  sea  level  or  at  altitude,  metabolic 
disturbances  in  acid-base  balance,  or  combinations  thereof).  It 
treats  the  chemical  buffering  and  gas  transport  systems  in  reason¬ 
able  detail,  including  both  Haldane  and  Bohr  effects,  and  it 
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reccgnizes  the  presence  of  many  transport  delays  that  are  themselves 
dependent  variables.  It  permits  convenient  exploration  of  a  variety 
of  possible  control  functions,  including  the  role  of  cerebrospinal 
fluid  (CSF)  hydrogen  ion  concentration,  and  of  C^-CC^  interaction  at 
the  peripheral  chemoceptors .  Some  of  these  features,  particularly 
the  buffer  equations  and  the  many  dependent  time  delays,  introduce 
formidable  computational  difficulties,  which  have  oe<~  overcome  in 
the  digital  simulation  herein  described. 


2.  GENERAL  DESCRIPTION  OF  THE  SYSTEM 

We  consider  two  major  components,  a  controlled  system  ('plant* 
or  'process'),  and  a  controlling  system  (.controller). 


2-1.  Controlled  System 

The  'plant'  comprises  three  major  compartments  (lung,  brain, 
tissue)  connected  by  the  circulating  blood  (Fig.  1),  The  brain 
compartment  is  separated  from  a  cerebrospinal  fluid  reservoir  by  a 
membrane  of  restricted  oermeability.  The  events  of  the  respiratory 
cycle  are  ignored  and  the  lungs  regarded  as  a  box  of  constant  volume, 
uniform  content,  and  zero  dead  space  ventilated  by  a  continuous  uni¬ 
directional  stream  of  gas*  If  the  alveolar  RQ  differs  from  unity, 
as  it  generally  will,  the  rates  of  flow  of  inspired  and  expired  gas 
will  differ.  We  assume  that  alveolar  gas  tensions,  PA(Cq  y  pA(q  y 
and  )  fire  equal  at  every  instant  to  those  in  expired  air  and 

arterial  blood  leaving  the  lung.  The  total  oxygen  content  of  this 
arterial  blood,  C^q  y  will  be  the  sum  of  physically  dissolved  and 
chemically  combined  (i-e.,  oxyhemoglobin)  components.  The  former 


varies  linearly  with  Po/n  ,,  whereas  the  latter  is  a  complex  function 

a 

of  both  )  ant*  arteri6^  blood  pH  (oxygen  dissociation  curve 

including  Bohr  effect  110J).  Similarly,  the  total  CC>2  content  of 


Fig.  1— - -The  controlled  system 
(Symbols  defined  in  App.  A.) 
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this  blood,  C  ,rn  ,  will  comprise  both  physically  dissolved  and 
chemically  combined  components,  the  former  varying  linearly  with 
Pa(co  )  and  t^ie  letter  being  a  complex  implicit  function  of  y 

the  standard  bicarbonate  content,  (BHCO^)^*  the  total  hemoglobin  and 


oxyhemoglobin  contents  (Hb,  HbC^),  and  the  plasma  protein  content 
(CC>2  absorption  curve  including  Haldane  effect  [10J).;  Carbamate  is 
not  explicitly  treated  as  a  separate  entity.  The  pH  of  this  blood 
will  depend  upon  its  ^  and  bicarbonate  content  as  defined  by 

the  Henderson-Hasselbalch  equation.  Only  physically  dissolved 
nitrogen  will  be  present. 

After  a  transport  delay,  which  depends  upon  vascular  volume  and 
blood  flow  rate,  this  arterial  blood  arrives  at  the  brain  or  tissue 
compartment.  We  regard  the  gas  tensions  in  each  of  these  reservoirs 
as  uniform  and  equal  to  those  in  exiting  venous  blood.  Oxygen  is 
removed  from  and  CO2  is  added  to  each  reservoir  at  a  constant  rate 


set  by  metabolism..  The  oxygen  and  nitrogen  contents  of  these  reser¬ 
voirs  are  taken  to  include  only  physically  dissolved  gas,  but  the 
total  CO2  content  will  include  a  bicarbonate  component  in  accordance 
with  the  CO2  absorption  curve  for  that  reservoir.  No  Haldane  effect 
is  attributed  to  the  reservoirs'  CC>2  buffers.  Although  we  assume 
reservoir-venous  blood  equality  for  CO2  tension,  the  bicarbonate 
and  hydrogen-ion  concentrations  will  not  necessarily  be  the  same. 

We  do  not  attempt  to  treat  this  "effective  impermeability"  in  any 
rigorous  physico-chemical  fashion  [11],  but  simply  reserve  the  right 
to  restrict  the  movement  of  certain  chemical  species  between  the 
reservoirs  and  the  blood..  In  this  way,  for  example,  we  can  simulate 
the  limited  availability  of  brain  and  CSF  bicarbonate  for  buffering 
stroig  a'.id  introduced  into  the  blood,  a  feature  of  the  "blood-brain 


barrier . 


It 


the  brain  compartment  communicates  with  a  cerebrospinal  fluid 
reservoir  via  a  membrane  permeable  to  respiratory  gases  only.  The 
latter  diffuse  across  this  membrane  at  rates  proportional  to  their 
tension  gradients.  The  entire  gradient  occurs  across  the  membrane, 
the  reservoir  concentrations  being  regarded  as  uniform.  The  CSF 
reservoir  differs  from  the  brain  comipartment  in  containing  no  protein 
capable  of  buffering  carbonic  acid.  For  our  purposes  we  regard  it 
simply  as  a  solution  of  sodicm  bicarbonate  whose  bicarbonate  content 
thus  remains  constant  at  all  CO2  tensions  above  10  mm  Hg  [10] . 

Venous  blood  leaving  the  brain  combines  with  that  leaving  the 
tissue  after  appropriate  time  delays  tc  form  mixed  venous  blood, 
which  after  another  delay,  enters  the  lung  to  complete  the  circle  of 
the  gas  transport  and  exchange  process.:  Total  cardiac  output  as  well 
as  the  local  blood  flow  to  the  brain  are  known  to  vary  as  functions 
of  arterial  CO2  and  oxygen  tensions  [12-16 ],;  Since  circulatory 
transport  delays  are  functions  of  blood  flow  rates,  these  delays 
become  dependent  variables  as  P  ,rn  \  and  P.  ,n  ..  change  in  response 
to  such  forcings  as  CO2  inhalation  or  hypoxia.  These  delays  are 
treated  as  pure  dead  times,  which  implies  a  square  blood  flow  front 
with  no  intravascular  "smearing.: " 


2.2.  Controlling  System 

This  includes  receptor  elements  that  monitor  concentrations  of 
certain  chemical  species  at  particular  locations,  afferent  nerves 
that  transmit  this  information  to  the  central  nervous  system,  the 
neural  centers  themselves,  motor  nerves  to  the  respiratory  muscles, 
the  muscles  u -...selves ,  and  finally  the  thorax-lung  pump  which  they 
drive.  As  in  our  original  model  as  well  as  those  that  have  followed 
it,  we  shall  not  attempt  to  describe  this  system  in  terms  of  the  details 
of  its  components  but  instead  will  go  direccly  from  chemical 


concentrations  at  receptor  sites  as  the  inputs  to  ventilation  as  the 
output..  In  so  doing,  we  shall  assign  no  dynamic  elements  to  this 
system,  thus  implying  that  any  delays  associated  with  it  are  neg¬ 
ligibly  jhort  compared  to  those  of  the  "plant."  This  is  certainly 
true  of  the  neural  components  and  probably  of  the  mechanical  elements 
as  well. 

3.  EQUATIONS  OF  THE  SYSTEM 

Using  symbols  and  units  defined  in  App.  A  and  "normal"  parameter 
values  given  In  App.  B,  we  first  write  material  balance  equations  for 
CO2,  C^,  and  N2  for  each  of  the  three  compartments,  lung,  brain,  and 
tissue: 


FA(C02) 

1 

[VIFI(C02)  VEFA(C02)  +  (#-47)  Q(Cv(C02)  Ca(C02)j 

]  (1.1) 

FA(02) 

1  j 

■  kl  1 

[VIFX(02)-VEFA(02)  +  (IS?)  0  (Cv  <  C2  )~Ca  <  C2  >)] 

(1-2) 

fa(n2) 

1 

"«L 

[Vl(N2)  VEFA(N2)  +  (b-47)Q(Cv(N2)  Ca(N2))] 

(1.3) 

CB(C02) 

-fc 

0 

[“^(COj)  +  QB(CaB(C02)~CvB(C02))_DC02(PB(C02)' 

’pcsf(co2))] 

r- 

(1-4) 

cb(02) 

1 

kb 

L_MRB ( 02 )  *  QB  VCaB(02)  CvB(02))~D02(PB(02)  PCSF(C>2))1 

d-5) 

CB(N2) 

1 

■«; 

_QB  (■ Ja  B  ( N2 ) ' ~CvB  ( N2  ))  ( PB  ( N2 )  ”PCS F  (N2 ))] 

(1-6) 

Sco2) 

1 

“ 

^®T(CU2)  +  (Q“V  (CaT(C02)~CvT(C02))] 

(1-7) 

CT(02) 

1_ 

Kt 

-MRi(o2)  +  (Q~Qb)  (CaT(02)“CvT(02))] 

(1-8) 

7- 


The  balance  (or  diffusion)  equations  for  the  CSF  subcompartment  are 
written  in  terms  of  gas  tensions: 


pcsf(co2)  ”  Cdco2  //KCSFka  CSF(C02))  CPB(C02)“PCSF(C02)) 
PCSF(02)  ■  CD02/KCSFkaCSF(02))  CPB(02)“PCSF(02)) 
PCSF(N2)  “  (DN2/KCSFkaCSF(N2))  CPB(N2)“PCSF(N2)) 


(1.10) 

(1.11) 

(1.12) 


By  adding  Fqs.  (1. 1)~(1. 3) ,  noting  that  the  volumetric  fractions 

of  C02,  02,  and  N2  must  add  to  one  in  both  inspired  and  alveolar  air, 

•  •  • 

and  that  the  sum,  FA^C0  ^  +  F^q  ^  +  FA^N  y  must  equal  zero  since 
the  lung  volume  is  constant,  we  can  express  V£  in  terms  of  Vjt 


VE  "  VI  +(f^7)  <?|(Cv(C02)~Ca(C02))  +CCv(02)“Ca(02))  +  CCv(N2)_Ca(N2) 

(2.1) 


We  now  define  alveolar-arterial  concentration  equilibria  for 
C02,  02,  and  N2*.  For  C02: 

Ca(CO.)  *  <BHC03>b  +  0-375[(Hb)-Cal  ^  J  -  [0.I6  2.3(Hb)] 


log 


'Ca(C02)-taC02<1M'7>F 


0.01(B-47)FA(CO  . 


a(co2)\ 

)_0-lj 


+  ,'aC05(B~47)FA(C0,) 


(3.1) 


For  02: 


Ca(0„)  "  ka0~(Br'47>FA(0  '  +  Ca(HbO~)J 


(3.2) 
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where 


Ca(Hb02)  '  <“»  L1  *  ^PC-S(B-47)FA(o2))]  , 

S  -  Q.44921(pH)  -  0.10098  (pH  )2  +  0.0066815  (pH  )3  -  0.454, 

g  a  a 


(3.3) 

(3.4) 


PHfl  -  9  -  log  Ca(H+), 

and 


(3.5) 


C 


a(H+) 


kaC02^B-47^FA(C02) 

-Ca  (C02)“kaC02^7>  FA  (C02  ) 


(3.6) 


Thus  Eq.  (3.1)  is  the  C02  buffer  equation  including  the  Haldane 
effect  [2,  10],  and  Eq.  (3.3)  is  the  02  dissociation  curve  including 
the  Bohr  effect.  Equations  (3.3)  and  (3.4)  were  fitted  empirically  to 
standard  data  on  human  blood  [17]. 

For  N2,  we  have  simply: 


Ca(N9)  "  kaN9(B^7)FA(N9) 


(3.7) 


We  next  define  venous  blood-brain  equilibria  for  C02,  02>  and  N2 
in  terms  of  gas  tensions.  For  C02: 


CB(C(V  "  (BHC°3>B  -  0  62 


/CB(C0,)  kaB(C02)PB(C0,) 
,108  V  O-01  PB(CO^ 


\  -  0.14 


*  kaB(C02) FB(C02)  ‘ 


(4.1) 


and 


CvB(C02)  ■  <BHC03>b  +  0  375  j<H1»-CvB(Hb02)]  -  [°'16  +  2'3<Hb>] 


los( 


CvB(CQ2)~kaC02PB(C02) 


1.01  P 


B(C02) 


+  kaC02PB(C02) ’  (4‘2> 


Equation  (4.1)  is  a  modified  C02  buffer  relation  for  brain  with 
the  Hb  and  Hb02  effects  removed,  and  Eq.  (4.2)  is  the  same  form  as 
Eq.  (3.1).  For  02: 


CvB(02)  "  aB(Q  }  CB(02)  +CvB(Hb02)’ 


(4-3) 


where  as  before 


CvB(Hb02)  "  <Hb>  [3  "  “',(-S(CB(02)/kog(02)>)]  - 


CB(02)  "  kaB(02)PB(02)» 


(4.4) 


S  -  0.44921(pHvB)  -  0.10098(pHvB)2  +  0. 0066815 (pHyg) 3  -  0.454,  (4.5) 


P«vB  "  9  ~  l08  CvB(H+) 


(4.6) 


CvB(H+ 


)  "  *'  b 


kaC02PB(C02) 
CvB(C02 )_kaC02PB ( C02 ) 


(4.7) 


Equations  similar  to  (4.6)  and  (4.7)  with  Cvg  replaced  by  Cg 
will  yield  Cg^-f^  and  pHg. 
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For  N2>  we  have  simply: 

“n2 

CvB(N2)  ’  aB(Hj)  CB(N2V 
CB(N2)  “  kaB(N2)PB(N2) 


(4.8) 


An  exactly  analogous  set  of  equations  with  values  for  the  tissue 
reservoir  (subscript  T)  and  tissue  venous  blood  (subscript  vT)  substi¬ 
tuted  for  those  of  the  brain  and  brain  venous  blood  is  used  to  define 
venous-tissue  equilibria  for  C02,  CU,  and  N2*  We  will  net  bother  to 
rewrite  them  but  simply  reserve  the  equation  numbers  [Eqs.  (5 . X) — (5 . 8) ] 
to  designate  them. 

We  note  here  that  the  hydrogen-ion  concentration  in  CSF,  which 
we  will  need  as  an  input  to  the  controller,  is  given  by: 


fkaCSF(C(v  )PCSF(G02)  1 
:CSF(H+)  “  K'[~  (BHC03)CSF  '  j 


(6.1) 


and 


pHCSF  "  9  ~  log  CCSF(H+). 


(6-2) 


We  next  define  the  dependence  of  cardiac  output  and  brain  blood 
flow  on  arterial  C C  and  02  tensions  using  information  from  the  lit- 
erature  [12-16]  and  assigning  an  arbitrary  first-order  lag  to  the 
responses.  For  cardiac  output,  we  have: 


Q  -  (fQ-  Q)/r,, 


(7-1) 


+  ^(co2)- 


where 


fQ  -  QN  +  AQ((j2) 


(7.2) 


11- 


* 

For  X  s  (B~47)F^q  )  *  ^a(0  V  anc*  *  <  we  ^ave; 

2  *  2 

AQ(02)  "  9-6651  -  0.2885X  +  (2. 9241E-3)X2  -  (1 . 0033E-5)X3  (7.3)* 

And  for  X  >  104: 

“)(o2)-°-  (?•*) 

For  Y  =  <B-47)Fa(COj)  -  pa(c02).  and  40  <  Y  <  60, 

aV>2)  -0.3(Y-40),  (7.5) 

And  for  fill  other  values  of  Y, 

6<5(C02)  ■  °-  <7-6> 

For  brain  blood  flow,  we  have: 

%  "  ^fQ  ~  QB^r2'  (7-7) 

where 

fQB  "  %N  +  AQB(02)  +  AQB(C02)’  <7,8) 

^"Floating  point"  notation,  i.e. :  E  -  n  s  x  lO"”*1,  in  this  and 
subsequent  equations. 


9 
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For  X  *  104: 

Qb(0  .  -  2.785  -  0.1323X  +  (2.6032E-3)X2 

-  (2. 324E-5)X3  +  (7 .6559E-8)}d ,  (7.9) 

And  for  X  >  104: 

4%(02)  ■  °-  <7'10> 

For  Y  <  38: 

AQB(C0  )  “  (2. 323E-2)  -  (3.1073E-2)Y 

+  (8. 0163E-4)Y2  (7.11) 

For  38  <  Y  <  44: 

AQb(co2)  -0,  (7-12) 

And  for  Y  >  44 : 

AQB(co  )  -  -15.58  +  0.7607Y  -  (1 . 2947E-2)Y2 

+  (9.3918E-5)Y3  -  (2. 1748E-7)Y4.  (7.13) 

In  each  case,  the  polynomials  in  X  or  Y  represent  empirical  fits  to 
data  taken  from  the  literature,  data  admittedly  incomplete  in  several 
details . 

We  now  express  the  arterial  concentrations  of  (X^,  O2,  and  N2 
at  the  entrance  of  the  brain,  and  of  the  tissue  reservoirs — in  terms 
of  their  concentrations  in  arterial  blood  leaving  the  lungs  at  an 
appropriate  time  in  the  past,  i.e.,  at  "lag  time"  (t  -  t),  where  t 
is  the  blood  transport  delay: 


V 
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CaB(C02)(t^  "  ^(CO^^aB^  (8’1* 
CaB(02)^  “  Ca(02)^t“TaB')  ^8'2^ 
CaB(N2)(t)  "  Ca(N2)^twraB^  (8-3> 
CaT(C02)^^  "  Ca(C02)^aT^  (8,4> 
CaT(02)(t)  "  Ca(02)(t^raT)  (8‘5) 
CaT(N2)(t)  "  Ca(N2)(t~TaT)  (8,6) 


and  in  similar  fashion,  we  express  the  mixed  venous  concentrations 
entering  the  lung  in  terms  of  concentrations  in  venous  blood  leaving 
the  brain  and  tissue  reservoirs  at  appropriate  lag  times: 

Cv(C02)(t)  "  ^[VCvB(C02)(t_TvB^  +  ^V(CvT(C02)^vT^] 

(8.7) 

Cv(02)(t^  "  Q[<VCvB(02)(t“TvB^  +  ^■%)(CvT(02)^C_TvT^] 

(8.8) 

Cv(N2)vt^  "  Q[QB^CvB(K2)^t-Tvb^  +  ((H^BHCVT(N2)^t~TvT^] 

(8.9) 

The  variables  defined  by  Eqs.  (8.1)-(8.9)  will  be  called  "lag 
variables"  for  convenience. 

Finally  we  note  that  the  transport  delays  are  not  constant  but 
vary  with  blood  flow  rates  as  defined  below: 
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1.062 


0.015 


aB 


t—r 


TaB~raB(l)  1  t-r 


aB(l) 


a  B 


-**  + 


(8.10) 


Q  dt 


■—*  f 


aB(l)  “  t—r 


aB(l) 


Qb  dt 


1.062 


0.735 


aT  r  ,t_raT(l) 


(8.11) 


aT  'aT(l)  “t-r 


I  Q  dt  -  j  (Q-Qg)  dt 

0  *. 1  /  1  \  <>  D 


aT 


aT(l)  t—r 


aT(l) 


0.06 


0.188 


Tvb 


1  r  VB(1)  ^  p 

— -  n_  Hr  — i — .  I 


(8.12) 


Qb  dt  J  Q  dt 


TvB~TvB(l/  J  t-irvB  ^BC1)  "^vBd) 


2.94 


0.188 


vT 


t-T 


i»TM  ■» 


(8.13) 


r-“ -  f  . <«->  dt  —4 -  j  Q  dt 

tvT^vT(1)  Jt-rvT  B  tvT(1)  Jt-rvT(1) 


1.062 


0.008 


ao 


Tao~Tao(l)  J t—r 


rH»(l) 

J  Q  dt 


(8.14) 


ao 


J^ao(I)  B 


Q*  dt 


The  last  value,  t„_,  which  does  not  appear  in  Eqs.  (8.1)-(8.9), 
is  the  lung-to-carotid  body  delay  and  will  be  used  in  the  controller 
function.  In  each  equation,  the  numerator  represents  the  volume  of 
the  appropriate  vascular  segment  through  which  blood  flows  at  the 
appropriate  ’'past-average  rate"  defined  by  the  denominator.. 
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We  have  now  written  all  of  the  relationships  required  to  define 
the  open-loop  operation  of  the  isolated  "plant. "  To  close  the  loops, 
it  remains  to  write  a  controller  function.  A  major  purpose  of  our 
simulation  is  to  explore  the  implications  of  a  variety  of  such  func¬ 
tions,  and  so  far  we  have  looked  at  two: 


VI  “  1,1  +  1-31  PB(C02)  +  TERM  VI(1)' 


(9-1) 


and 


VI  *  1,138  CCSF(H+)  +  1,154  ^(H+^^ao^  +  TERM  VI(2) 

l4.9 
JV 


TERM  =  (23.6E-9)[l04-(B-47)FA(0  )(t-rao) 


for  (^7)FA(02)(t-rao)  <  104, 


TERM  s  0 


for  (B^7)FA(02)(t-rao)  >  104, 


(9.2) 


(9-3) 

(9-4) 


and  are  constants  so  adjusted  that  P^^)  at  rest» 

breathing  air  at  sea  level  (standard  control  point).  Equation  (9.1) 
is  similar  to  that  used  in  our  earlier  model,  whereas  Eq.  (9.2) 
incorporates  recent  findings  on  the  roles  of  CSF  and  peripheral 
receptors  [18-21].  Neither  includes  an  02-C02  interaction  term,  but 
this  as  well  as  a  variety  of  other  possibilities  can  easily  be 
introduced. 


4.  SOLUTION  OF  SYSTEM  EQUATIONS  BY  DIGITAL  COMPUTER 

The  equations  were  programmed  for  solution  on  IBM  /0 44  and 
CDC  3400  digital  computers.  For  this  purpose,  the  basic  material 
balance  equations  (Eqs.  (1.1 )— ( 1 • 9 ) ]  were  rewritten  as  differential- 
difference  equations  with  variable  time  delays,  using  the  appropriate 


interrelationship  equations  developed  above: 


a<co2)  -  vifi(C02)  vefa(C02>  4-  (^7)(qbcvB(Co2) 


+  (Q“QB)CvT^C0  ^(t-rvT)  Q  Ca(co  y 


(t^vB> 


(10.1) 


1  ..  r-  r.  ,  '*63 


a(°2)  "  IFI(02)  VEFA(02)  +  lb-47AQBCvB(02)(twrvB) 

+  (^-QB)CvT(02^(t-rvT)  -  Q(kao2(Er^7)FA(02)  +  Ca(Hb02))) 


A(N9) 


K^I*!^)  VEFA(N2)  +  (^47)  (QBCvB(N2)(t"-rvB) 


(10.2) 


+  <^VCvT(N2)(t^vT>  "  QK  (B-47)pA(N  )))]  (10.3) 


'B(C0 


;)  "  ^[^(COg)  +  QB(Ca(C02) 


(t^TaB)  CvB(C02) 


“cojC1 


'2VPB(C02)  PCSF(C02)/ 


(10.4) 


,)  "  ^[^(O^  +  QB^ka02(B,~47)FA(02)  +  Ca(Hb02))  (t_TaB) 

~  (a^y)CB(02)~CvB(Hb02))  “  D02(PB(02)_PCSF(02))j 
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:B(N2)  *  Kg  QB\JcaN2(Bf*47^FAvN2)(t:  '  aB* 


(^t-) 


JB(N,) 


-  %(' 


vPB(N2)  PCSF(N2) 


(10.6) 


JT  ( CO 


02)  *  K.r(^®T(C02)  +  (Q_<V  (C 


^(CO^^aB*  CvT(C02) 


(10,7) 


jt(o2) 


Kt[^  MRt(°2)  +  (Q~QB^{(ku02(Br^7)FA(02)  +  Ca(Hb02))(t  TaT^ 

■  (^fr)  CT(02rCvT(Hb02)  }  (1° 


(10.8) 


jT(N9) 


((Q-Qb)/^)  K2(B^7)Fa(N2))  (t-raT) 


/  QN„  \ 


Ct(N2>J 


(10.9) 


Equation  (2.1)  defining  Vg  in  terms  of  Vj  was  also  rewritten  to 
incorporate  the  appropriate  lag  times: 


+  (f!t)  \^BL(CvB(ro2)  +  CvB(02)  +  CvB(N2))  ^"^vB^] 
(Q_QB)[(CvT(C02)  +  (i^))CT(02)  +  CvT(HbC2)  + 

Q[Ca(C02)  +  ka02(Br  ;)FA(02)  +  Ca(Hb02)  +  kaN2(B_47 '*  FA(N2)]) 


(11,1) 


0 


-i&~ 


Methods  for  the  numerical  integration  of  such  differential- 
difference  equations  with  several  time  lags  have  been  considered  in 
detail  elsewhere  [22-24].  Here  we  shall  only  point  out  certain 
special  features  of  the  present  program  that  deserve  comment. 

The  first  concerns  the  calculation  of  current  values  for 

Ca{C02)’  Ca(Hb02)’  CvB(C02)‘  CvB(Hb02)»  CvT(C02)’  and  CvT(Hb02)  t£> 
put  into  this  system  of  equations.  Special  problems  are  involved 

here  because  of  the  interdependence  of  CO 2  content,  oxyhemoglobin 

content,  and  pH,  and  because  of  the  fact  that  the  buffer  equation 

cannot  be  solved  explicitly  for  the  desired  variable.  Let  us  examine 

the  procedure  used  here.  Starting  with  the  current  value  of  F./r>r.  v 

A(C02) 

and  a  trial  guess  for  Cfl^co  ^  (say  Ca^CQ  we  calculate 

Ca(H+)(l)  and  pHa(l)  ^roin  Ecls'  (3*6)  and  (3.5).  Putting  this  value 
of  into  Eq.  (3.4)  and  using  the  current  value  of  F^q  y  we 

obtain  from  Eq.  (3.3).  Then  using  y  we  solve 

Eq.  (3.1)  for  )(2)‘  by  a  Process  of  iterative 

guessing,  putting  ca(Co2)(l)  on  the  ri8ht  hand  side,  calculating 
Ca(CO  ) (2)  on  the  leftj  comparing  the  two,  adjusting  Ca^0  and 

repeating  until  the  two  agree  within  a  predetermined  limit  of  error. 
Now,  using  this  new  value  for  Cfl^Q  y  we  repeat  the  entire  process, 
recalculating  pHa,  and  C,^  until  we  get  a 

that  satisfies  Eq.  (3.1)  on  the  first  trial.  Going  next  to  the  brain, 
we  solve  Eq..  (4.1)  for  P^CC^)’  us*n§  t*le  current  value  of 
Again  we  use  an  iterative  guessing  procedure,  but  since  Hb02  is  not 
involved  here,  it  is  unnecessary  to  include  the  pH  and  Hb02  steps.. 

Then  using  this  value  of  PB(C0  the  current  value  of  Cg^0  y  and 

an  initial  trial  guess  for  )>  we  apply  the  complete  iterative 

procedure  to  Eqs <  (4.2)-(4,7)  to  obtain  current  values  for  CvB(C0  y 
CvB(0  )  and  p51vB'  FinaHy>  going  to  the  tissue  reservoir,  we  apply 
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the  same  technique  to  Eqs .  (5.1)— (5.7)  to  obtain  current  values  for 

PT(C02r  CvT(C02)’  CvT(02)’  and  pHvT' 

A  second  feature  involves  the  storage  and  retrieval  of  the  "lag 

variables"  defined  by  Eqs.  (8.1)-(d.9).  Provision  is  made  for  storing 
in  the  computer  previous  values  of  these  variables  in  a  table,  con¬ 
tinuously  updated,  that  extends  an  appropriate  distance  into  the  past, 
usually  0.5-1. 5  minutes.  After  calculating  current  values  for  the 
lag  times,  an  interpolation  routine  is  used  to  retrieve  the  proper 
values  of  the  lag  variables  for  use  in  the  next  integration  step. 

In  early  versions  of  the  program,  current  rather  than  average  blood 
flow  values  were  used  to  calculate  the  required  t  values,  but  in 
later  versions,  the  correct  average  values  were  computed.  To  obtain 
these  it  was  necessary  to  integrate  the  various  blood  flow-time  func¬ 
tions  "backwards"  into  the  past,  starting  at  current  time  and  stopping 
when  the  value  of  the  integral  became  equal  to  the  appropriate  vascular 
volume.  This  point  defined  the  lower  limit  of  the  integration  and 
hence  the  appropriate  t  . 

Finally,  for  convenience  in  future  explorations  and  modifications, 
the  programs  have  been  written  in  the  form  of  multiple  subroutines, 
each  of  which  can  be  altered  more  or  less  independently  of  the  others. 
The  latest  version  of  our  Fortran  IV  program  includes  approximately 
500  statements.*  A  sample  tabular  output  is  provided  in  App.  C,  and 
graphs  of  the  material  will  be  discussed  below. 

5 .  RESULTS 

The  equations  of  the  model  and  the  corresponding  digital  simula¬ 
tions  have  evolved  through  several  stages  cf  increasing  complexity. 

In  the  first  version  there  was  no  CSF  compartment,  an'4  values  were 

Programs  are  available  from  the  authors  on  request. 


0 


computed  from  current  rather  than  "past-average”  blood  flows.  The 

A 

controller  Inputs  were  brain  P ^  ,  brain  (H  ),  and  arterial  PQ  at 

2  2 

the  peripheral  chemoceptors .  A  variety  of  forcings  were  studied 
including:  (a)  COj  inhalation  (1,  3,  5  percent)  and  recovery;  (b) 
hypoxia  at  sea  level  (10,  8,  6  percent  O2)  and  recovery;  and  (c) 
altitude  hypoxia  (15,000,  18,000,  20,000  feet).  With  few  exceptions, 
the  solutions  were  stable,  and  the  quasi-steady  state  values  were 
close  to  those  expected  in  the  prototype  (Table  1).  It  is  noteworthy 
that  a  true  steady  state,  if  judged  by  the  alveolar  RQ  calculated 
at  each  step,  was  not  attained  in  30  minutes  for  any  of  the  forcings 
except  perhaps  1  percent  CC^.  Transient  behavior  during  CO2  in¬ 
halation  and  recovery  was  similar  in  form  to  that  of  an  earlier 
model  [3];  an  example  is  given  in  Fig.  2.  However,  it  will  be  noted 
that  the  "on  response"  of  ventilation  is  much  faster  in  the  present 
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Fig.  2 — Response  of  Version  1  to  5%  CO2  "pulse.” 
Pulse  on  at  zero  time,  off  at  30  min. 
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model,  requiring  only  3  minutes  to  reach  the  80  percent  level  rather 
than  the  12  minutes  needed  by  the  original  version  [3].,  This  dif¬ 
ference  has  also  been  noted  by  others  (5,  9]  who  have  separated  out 
a  brain  compartment  from  the  single  lumped  tissue  reservoir  previously 
assumed  [3];  it  is  directly  related  to  the  relative  ratios  of  com¬ 
partment  volume  to  blood  flow.  Closer  attention  to  the  early  details 
of  the  prototype  response  [19-21]  has  indicated  that  neither  model  is 
too  faithful,  for  there  appears  to  be  an  initial  fast  response  of 
relatively  modest  magnitude  followed  by  a  much  slower  but  larger 
component.  Moreover,  the  ventilatory  response  lags  considerably 
behind  the  change  in  jugular  venous  after  the  first  minute  [19]. 

These  matters  will  be  referred  to  again  later.  As  expected,  a  short 
interval  of  periodic  breathing  occurred  on  sudden  return  to  air  after 
a  prolonged  period  of  hypoxia  (Fig.  3).  Figure  4  shows  a  phenomenon 
occasionally  j-ncountered,  i.e.,  the  sudden  appearance  of  numerical 
instability  in  a  previously  stable  solution  in  this  case  after  35 
minutes  at  an  altiti  de  of  20,000  feet.  Disappearance  of  this  behavior 
on  changing  the  integration  step  size  confirms  the  suspicion  that  it 
is  a  computational  and  not  a  model  instability.* 

In.  the  second  stage  of  development,  a  CSF  reservoir  was  added, 
but  t  values  were  still  computed  from  current  blood  flows.  The 

X  X 

controller  inputs  were  now  the  CSF  (H  ) ,  and  both  arterial  (H  )  and 

With  single  precision  arithmetic,  a  step  size  of  0.01  minute 
was  satisfactory  on  the  CDC  3400  (48-bit  word).  On  the  IBM  7044 
(36~bit  word),  a  step  size  of  0.007812,  i.e.,  1/27,  was  usually 

Q 

satisfactory,  but  occasionally  a  smaller  one  was  required  (1/2 
or  1/29). 


TIME  -  MINUTES 

Fig.  3 — Periodic  breathing  on  recovery  from  hypoxia  (30  min  on  10%  0 
Version  1).  Value  )  changed  exponentially  from  0.10  to  0.21 

‘during  first  min. 


TIME  -  MINUTES 

Fig.  4 — Sudden  appearance  of  numerical  Instability  after  35  min  at 

20.000  ft.  Version  1. 
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Pn  at  the  peripheral  chemoceptors,  a  control  scheme  similar  to  that 
°2 

suggested  by  Lambertsen  [19],  and  compatible  with  the  recent  work  of 
MLtehell,  et  al.  [183-  This  control  function  should  introduce  both 
fast  and  slow  components  into  the  CO2  response,  and  should  also  permit 
a  preliminary  study  of  responses  to  metabolic  acidosis.  Figure  5 
shows  the  behavior  of  ventilation,  CSF  (H  ),  and  arterial  (lung  exit) 
(H+)  during  and  following  a  30  minute  "pulse"  of  5  percent  COj  in 
inspired  air.  The  similarity  of  the  "on  transient"  to  the  behavior 
observed  in  dogs  by  Lambertsen  (19]  is  apparent.  In  contrast  with 
Version  1,  ventilation  once  more  requires  12  minutes  to  reach  80 
percent  of  its  total  response,  but  in  contrast  with  our  original 


Fig.  2 — Response  of  Version  1  to  5%  CO2  “pulse." 
?ulse  on  at  sero  time,  off  at  30  min. 
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model  [3]  and  in  agreement  with  recent  observations  [19—21],  we  now 
also  have  an  initial  rapid  component  so  that  about  30  percent  of 
the  response  occurs  within  one  minute.  The  exact  relationship 
between  the  three  curves  in  Fig.  5  can  be  altered  by  changing  the 
relative  weightings  assigned  to  peripheral  and  CSF  (H+j  receptors. 
Although  the  entire  time  course  of  the  response  to  metabolic  acidosis 
was  obtained,  the  transient  portion  is  not  really  very  relevant,  for 
the  forcing  was  introduced  in  a  highly  artificial  manner,  i.e., 
a  step  change  in  (BKCO^)^  and  (BKCO^)^  at  zero  time  while  keeping 
(BHCO^g  and  (BHCO^^gp  constant.  However,  it  is  of  interest  that 
a  stable  solution  was  obtained  and  that  the  steady  state  values 
shown  in  Table  2  are  quite  realistic  [2].  Thus,  the  same  reduction 
in  arterial  pH  produced  by  5  percent  C09  inhalation  on  the  one  hand, 
and  by  metabolic  acidosis  on  the  other  was  associated  with  a  ventila¬ 
tion  ratio  of  3.66  in  the  former  and  only  1.25  in  the  latter.  Whereas 

—  .  ...  .  ..  .  — 4*. 

lu 2  innaiacion  proaucea  nypercapxua  ana  mereasea  (h  )  cnrougnouc 

a- 

the  body,  the  increase  in  (H  )  in  metabolic  acidosis  was  confined 
to  the  blood  and  "not-brain"  tissue,  while  the  (H+)  of  brain  and 
CSF  actually  fell.  Thus  the  general  qualitative  and  quantitative 
behavior  of  the  model  seems  quite  reasonable  under  these  conditions. 

The  major  modification  introduced  into  the  third  version  of  the 
model  was  the  use  of  correct  "past-average"  blood  flow  values  to 
compute  current  values  for  the  time  lags.  This  did  not  produce  any 
significant  difference  in  the  response  to  a  5  percent  CO2  pulse. 

The  computer  generated  graph  of  Fig.  6  includes  the  time  course  of 
brain  (or  internal  jugular)  (H+)  which  is  identical  with  that  of 
brain  Pni~  on  this  non-dimensional  plot,  and  this  leads  the  ventila— 

Wn 

jC 

tory  response  after  the  first  minute  of  the  "on"  transient,  in 
agreement  with  experimental  observations  noted  previously  [25]. 
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Table  2 

STEADY  STATE  VALUES— VERSION  2 


Variable* 

Normal 

Control 

Metabolic 

Acidosis 

5% 

co2 

VI 

5.37 

6.71 

19.66 

VE 

5.33 

6.67 

19.62 

FA(C02) 

.0526 

.0421 

.0644 

fa(02) 

.1514 

.1632 

.1941 

FA(N2) 

.7959 

.7947 

,7415 

PA(C02) 

37.5 

30.0 

45.9 

PB(C02) 

47.8 

43.4 

55.3 

PCSF(C02) 

47.8 

43.4 

55.2 

PT(C02) 

42.3 

34.5 

50.1 

pa(02) 

108.0 

116.3 

138.4 

pb(02) 

36.5 

32.4 

45.0 

PCSF(02) 

36.2 

32.1 

44.6 

pt(02) 

46.7 

50.6 

57.0 

Ca(C02) 

.5629 

.3966 

.6048 

CB(C02) 

.6396 

.6185 

.6718 

.6172 

.6142 

.6222 

CSF(COo) 

CvB(C02) 

.6309 

.4855 

.6588 

CT(C02) 

.6130 

.4212 

.6500 

CvT(C02) 

.5975 

.4300 

.6314 

Ca(Hb02) 

.1981 

.1982 

.1994 

CvB(Hb02) 

.1324 

■  1120 

.1483 

CvT(Hb02) 

.1592 

.1608 

.  1705 

PHa 

7.427 

7.370 

7. 368 

PHB 

7.375 

7.404 

7 . 331 

1 • 


* 


» 
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Ear  ly  details  of  the  "on"  and  "off"  transients  appear  in  the  expanded 
time  scale  graphs  of  Fig.  7.  The  response  to  a  30  minute  "pulse"  of 
10  percent  oxygen  with  subsequent  return  to  air  is  shown  in  Figs.  8-10. 
Prominent  features  are  the  ventilation  overshoot  during  the  "on" 
transient  (Fig.  8),  and  the  prolonged  interval  of  periodic  breathing 
during  recovery  (Fig.  10).  This  periodicity  persists  longer  than  it 
did  in  Version  1  (Fig.  3),  but  in  the  latter,  return  to  air  was  accom¬ 
plished  gradually  over  one  minute  whereas  in  the  present  example 
it  was  instantaneous.  Also  nicely  illustrated  is  the  sequence: 

PA(0  ^  fall— delay-ventilation  rise-PA^Q  ^  and  fall  at  the 

start  of  hypoxia  (Fig.  8),  as  well  as  the  sequence:  rise- 

delay— ventilation  fall-PA^CQ  ^  and  Ca^H+j  rise  at  the  start  of 
recovery  (Fig.  9). 


6 •  DISCUSSION 

A  major  purpose  of  the  present  study  was  to  provide  a  suf¬ 
ficiently  detailed  mathematical  description  and  digital  simulation 
of  the  "plant"  so  that  a  significantly  wide  variety  of  forcings 
and  control  hypotheses  could  be  explored.  This  has  been  accomplished 
to  some  extent,  the  plant  description  taking  the  form  of  a  set  of 
differential-difference  equations,  incorporating  a  number  of  dependent 
time-delays,  which  express  the  basic  material  balance  relationships 
of  the  system.  Additional  equations  are  required  to  define  the 
chemical  details  of  gas  transport  and  acid-base  buffering  mechanisms, 
concentration  equilibria  in  various  parts  of  the  system,  and  the 
dependence  of  blood  flow  rates  upon  CO2  and  Oj  concentrations. 

Finally,  a  control  loop  manipulating  pulmonary  ventilation  according 
to  any  one  of  n  variety  of  possible  schemes  can  be  included.  The 
model  is  sufficiently  general  to  accept  many  sorts  of  forcings,  e.g., 
CO2  inhalation,  hypoxia  at  sea  level  or  at  altitude,  metabolic 
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Fig.  9— Computer  plot  of  behavior  of.  Version  3 
during  first  0.75  min  following  instantaneous  return  to 
air  after  30  min  on  10X 


Fig,  10 — Continuation  of  Fig.  9  with  compressed  time  scale  to  show 
subsequent  course  of  recovery  from  hypoxia. 
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acidosis  and  alkalosis,  or  any  combinations  thereof,  and  any  arbitrary 
time  course  can  be  assigned  to  them.  Concentrations  of  a  variety 
of  chemical  species  are  available  at  several  different  locations 
for  use  in  possible  control  functions.  Extension  of  the  model  to 
include  any  number  of  additional  distinct  tissue  compartments  is 
simple  in  principle,  as  in  the  inclusion  of  other  inhaled  agents 
such  as  anesthetics. 

Solution  of  this  system  of  nonlinear  equations  would  be  completely 
impractical  without  computer  assistance.  The  digital  computer  is 
particularly  suited  for  large  multivariate  systems  incorporating 
time  delays,  although  the  storage  requirements  for  retaining  the  past 
histories  of  many  variables  may  get  uncomfortably  large.  To  our 
knowledge,  this  is  the  first  time  that  a  set  of  differential— difference 
equations  whose  time  delays  themselves  are  dependent  variables  has 
been  solved  by  digital  techniques  [22]-  The  latest  version  of  cur 
program  requires  about  one  minute  of  machine  time  (CDC  3400)  for 
every  four  minutes  of  experimental  simulation. 

Indicative  of  the  importance  of  dynamic  simulation  is  the  fact 
that  many  recent  contributions  to  our  understanding  of  the  "chemical 
control"  of  ventilation  have  been  based  upon  the  study  of  dynamic 
transients  in  an  attempt  to  separate  peripheral  (fast)  from  central 
(slow)  components  [19—21].  Such  studies  led  Lambertsen  [19]  to 
suggest  a  dual  sensing  mechanism  for  the  ventilatory  response  to 
C02»  one  receptor  responding  directly  to  arterial  (H+)  and  the  other 
to  the  (H  )  of  an  appropriate  central  fluid,  perhaps  CSF.  This 
suggestion  fits  in  very  well  with  the  results  of  Dejours'  "two  bieath 
CO?  test"  before  and  after  chemoceptor  denervation  [20,  21]  as  well 
as  with  other  observations  indicating  the  existence  of  superficial 
medullary  receptors  sensitive  to  CSF  (H  )  [18].  A  control  scheme 
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of  this  sort  incorporated  into  our  model  closely  reproduced 
Lambertsen's  observations.  In  both  experiment  and  model  (Fig.-  6), 
arterial  PCQ  (or  (H*))  changed  most  rapidly,  followed  in  order  by 
brain  venous  F^g  (or  (H  )),  ventilation,  and  finally  CSF  Pco 
(or  vH  )).  It  is  not  essential  to  this  scheme  that  the  central 
receptors  be  superficial  and  exposed  directly  to  "bulk  CSF,"  for 
they  could  lie  somewhat  deeper  within  the  medulla  in  accordance 
with  the  classical  view  and  the  recent  studies  of  Pappenheimer, 
et  al.  [25],  responding  to  the  (H+)  of  the  appropriate  local  central 
fluid.  It  is  essential,  however,  that  there  be  a  separate  arterial 
receptor  for  (H+)  else  the  results  of  chemoceptor  denervation  could 
not  be  accounted  for  [20,  21}. 

Dejours  [20],  and  Bouverot,  et  al.  [21],  have  also  used  a 

"two  breath  C>2  test"  to  study  the  peripheral  (fast)  response  to  100 

percent  ©2*  They  noted  the  difficulty  of  assessing  c'ne  mechanism  of 

the  ventilatory  response  to  long  lasting  pure  oxvgen  inhalation 

because  of  the  operation  of  many  other  factors  involving  the  gas 

transport  systems  Such  factors  are  taken  into  account  in  the 

present  model,  and  prolonged  administration  of  100  percent  O2  in 

Version  1  resulted  in  a  mild  hyperventilation  associated  with 

arterial  hypocapnia  and  alkalemia  together  with  brain  and  venous 

hypercapnia  and  increased  (H+)  [26].  The  initial  depression  of 

ventilation  observed  by  Dejours  was  not  reproduced,  for  we  have 

not  yet  rewritten  our  control  function  to  include  an  "O2  threshold" 

as  high  as  200  ran  Hg.  We  are  confident  that  this  modification 

together  with  an  appropriate  interaction  term  involving  arterial 

(H+)  and  P0  at  the  carotid  chemoceptors  will  permit  Version  3  of 
°2 

our  model  to  reproduce  Dejours'  CO2  and  09  results  in  detail.; 
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Although  step  or  pulse  inputs  of  gases  are  convenient  and  widely 
used  forcings,  che  production  of  metabolic  disturbances  in  acid-base 
balance  is  carried  out  in  a  variety  of  nonstandard  ways  ranging  from 
intravenous  injection  of  acid  or  alkali  to  prolonged  ingestion  of 
NH^Cl  or  NaHCO^  or  diabetic  sccuraula tit n  of  ketone  bodies.-  The  nature 
of  the  transient  response  is  of  course  critically  dependent  upon  the 
time  course  of  the  forcing  function,  and  while  it  is  possible  to 
simulate  a  wide  variety  of  these  in  the  model,  we  have  not  yet  carried 
out  any  such  systematic  exploration.  Very  recently,  the  primary  role 
of  the  peripheral  cherooceptors  in  the  ventilatory  responses  to  metabolic 
acid-base  disturbances  as  conceived  by  Mitchell,  et  al.  [27],  has 
been  questioned  by  Fend,  Miller,  and  Pappenheimer  [28],  who  found 
changes  in  CSF  (H+)  deemed  sufficient  to  account  for  the  behavior  of 
ventilation.  There  also  seems  to  be  disagreement  as  to  wheti er 
chemoceptor  denervation  abolishes  ventilatory  responses  to  metabolic 
acid— base  disturbances  [29,  30].  Although  some  of  this  disagreement 
may  be  more  apparent  than  real  and  dependent  upon  such  factors  as 
the  severity  and  duration  of  the  disturbance,  pre^e^te  or  absence  of 
anesthesia,  species  difference,  etc.,  it  is  apparent  that  much  more 
work  is  needed  here  to  define  the  details  of  material  transport  and 
buffering  mechanisms  in  the  brain  as  well  as  the  location  of  the 
central  chemoceptors . 

Despite  some  obvious  defects,  the  present  model  thus  seems  to 
provide  a  useful  dynamic  summary  and  synthesis  of  much  of  our  present 
understanding  of  lung-blood-tissue  gas  transport,  acid-base  buffering, 
and  what  is  generally  called  the  "chemical  control  of  pulmonary 
ventilation."  Although  many  of  its  features  can  also  be  found  in 
other  recent  models  [5-9],  the  present  version  is  perhaps  somewhat 
more  complete  and  general.  It  can  serve  as  a  convenient  base  for 


studying  sensitivity  to  parameter  variations,  for  exploring  modified 
control  modes,  and  for  gradual  expansion  to  include  aspects  so  far 
neglected.  With  this  in  mind,  the  digital  simulation  program  has 
been  made  flexible  enough  to  allow  convenient  modification  and 
expansion. 

That  both  modification  and  expansion  will  be  necessary  is 
abundantly  clear:  <Xir  treatment  of  tissue  buffering  and  material 
exchange  across  the  blood-brain  barrier  is  at  best  only  a  very  crude 
first  approximation.  As  noted  above,  detailed  knowledge  of  these 
matters  is  both  important  and  lacking.  The  model  does  not  account 
for  the  hyperpnea  of  exercise,  and  although  we  could  easily  introduce 
a  term  proportional  to  ©2  consumption  into  our  control  function  to 
accomplish  this  [31],  the  "physiological  isomorphism"  of  tach  a  term 
still  seems  to  be  obscure  [32].  Moreover,  attempting  to  account  for 
exercise  hyperpnea  in  this  way  would  make  it  painfully  obvious  that 
the  model  Includes  no  really  general  description  of  the  control  of 
cardiac  output  and  regional  blood  flow.  Since  the  respiratory  and 
cardiovascular  systems  are  really  part  of  a  single  gas  transport 
system,  close  coupling  between  their  control  mechanisms  might  be 
expected,  and  it  would  seem  that  a  real  understanding  of  either 
requires  an  understanding  of  both  [33]. 

Coupling  between  various  physiological  systems  treated  as  a 
collection  of  nonlinear  oscillators  has  been  considered  by  Iberall 
and  Cardon  [34],  and  spectral  analyses  of  skin  temperature  [34]  and 
ventilation  [34,  35]  time  series  obtained  from  normal  resting  men 
have  been  carried  out  in  efforts  to  identify  significant  frequencies. 
Four  major  periods  appear,  which  Iberall  suggests  may  be  associated 
with  a  metabolic  rate  cycle  (1-2  minutes),  a  vasomotor  cycle  (5~10 
minutes),  a  tissue  gas  storage  cycle  (20-40  minutes),  and  a  tissue 
heat  storage  cycle  (2  hours).  Our  model  includes  tissue  gas  storages 


-39- 


with  relatively  long  time  constants,  short  legs  in  the  response  of 
cardiac  output  and  brair  blood  flow,  and  a  number  of  short  blood 
transport  delays.  The  latter  are  primarily  responsible  for  the  short 
period  oscillations  (about  0.5  minute)  observed  in  the  model  during 
the  early  stages  of  recovery  from  hypoxia  or  CO2  inhalation.  It  is 
certainly  possible  that  tissue  metabolic  rate  is  normally  periodic 
at  rest  as  suggested  by  Iberall,  and  it  is  well  known  that  capillary 
blood  flow,  at  least  in  some  vascular  beds,  is  intermittent  [36], 
but  we  have  not  yet  included  these  possibilities  in  our  model. 

Finally,  although  the  "chemorespiratory"  physiologist  interested 
primarily  in  the  physics,  chemistry,  and  control  of  lung-blood— tissue 
gas  transport,  buffering  and  exchange  might  be  reasonably  happy  with 
this  model,  the  "neurorespiratory"  physiologist  interested  primarily 
in  the  intimate  central  neural  mechanism  of  respiratory  cycle  genera¬ 
tion  will  not  be  happy  with  it  at  all.  Certainly  there  are  no  illumina¬ 
ting  insights  here  for  him!  Perhaps  this  points  out  the  most  fruitful 
direction  for  further  worl  both  experimental  and  theoretical,  providing 
we  remember  that  ’’neural  control"  and  "chemical  control"  are  not 
mutually  exclusive  but  are  interrelated  parts  of  a  single  mechanism. 
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Appendix  A 

SYMBOLS  AND  'INI  ’S  FOR  EQUATIONS  OF  SYSTEM 


A.l.  EXPLANATION 

A  complete  symbol  has  to  general  form.  j,  where  A  designates 
the  class  of  variable  (e.g.,  solubility  coefficient,  concentration, 
metabolic  rate,  etc.),  i  spe  ifies  the  location  (e.g.,  brain,  blood 
at  lung  exit,  etc.),  and  j  i  ^ntifies  the  chemical  species  (e.g., 

-f- 

O2,  CO2,  H  ,  etc.).  In  cert  in  cases,  one  or  both  subscripts  are 
omitted  (e.g.,  Q  s  cardiac  c  itput;  B  =  barometric  pressure). 


A. 2.  TABLE  OF  SYMBOLS  AND  UNITS 


Subscripts 

Char¬ 

acter 

Loca¬ 

tion 

(i) 

Species 

(j) 

Definition 

Units 

a 

solubility  coefficient 

•  gas 

for  gas  in  blood 

liter  (STPD) /liter 

blood/atm,  37°C 

B 

(gas) 

for  jas  in  bn  in 

liter  (STPD) /liter 

brain/atm,  37°C 

T 

(gas) 

for  gas  in  tissue 

liters  (STPD) /liter 

tissue/atm,  37°C 

CSF 

(gas) 

for  ’.as  in  cerebro- 
sp  nal  fluid 

liters  (STPD) /liter 
CSF/atm,  37°C 

B 

barometric  pressure 

nm  Hg 

(BHCO3) 

stenda  d  bicarbonate 
content 

b 

of  Hood 

liters  C02(STPD)/ 

liter  blood,  37°C 

B 

of  brain 

liters  C02(STPD)/ 

liter  brain,  37 °C 

T 

of  tissue 

liters  C02(STPD)/ 

liter  tissue,  37°C 

0 
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Subs 

crlpts 

f 

i 

Char¬ 

acter 

■CTrta 

1 

Species 

U) 

Definition 

Units 

CSF 

of  cerebrospinal  fluid 

liters  C02(STPD)/ 

liter  CSF,  3?°C 

C 

concentration 

(ga°) 

of  gas 

liters  (STPD) 

(HbO.?) 

(H+) 

of  oxyhemoglobin 
of  hydrogen  ion 

liters  02(STPD) 
nanomoles 

a 

in  blood  at  lung  exit 

liter  blood-* 

aB 

in  blood  at  brain 
entrance 

liter  blood  * 

aT 

in  blood  at  tissue 
entrance 

liter  blood  * 

B 

I 

in  brain 

liter  brain-1 

CSF 

, 

in  cerebrospinal  fluid 

liter  CSF-1 

T 

in  tissue 

liter  tissue-1 

V 

in  blood  at  lung 
entrance 

liter  blood  1 

ao 

in  blood  at  carotid 
body 

liter  blood  1 

vB 

in  blood  at  brain 
exit 

liter  blood  1 

vT 

in  blood  at  tissue 
exit 

liter  blood  1 

D 

gas 

diffusion  coefficient 
for  gas  across 
"blood-bra in"  barrier 

liters  (STPD) /min/ 
mm  Hg 

F 

volumetric  fraction 

dimensionless 

(gas) 

of  gas 

A 

in  dry  alveolar  gas 

E 

in  dry  expired  gas 

I 

in  dr^,  inspired  gas 

{Hh 

blood  oxygen  capacity 

liters  (STPD) /liter 
blood 

"V. 

! 

conversion  factor 

atm/mm  Hg 

K’ 

I 

dissociation  constant 

i  1 

for  carbonic  acid 

na nomoles/liter 

K 

i 

i 

; 

lUrae 

liters 

O 

o  b^ain 

1 

! 

cs«- 

1 

' 

.  c-  ri  jospinal  fluid 
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Subscripts  1 

Char¬ 

acter 

HU  ■ 

Species 

<J> 

Definition 

Units 

T 

of  tissue 

L 

of  alveoli 

MR 

metabolic  rate 

liters  (STPD)/min 

(co2) 

of  carbon  dioxide 
production 

(C^) 

of  oxygen  consumption 

B 

by  brain 

T 

by  tissue 

P 

tension  or  partial 
pressure 

mm  Hg 

(gas) 

of  gas 

A 

in  alveoli 

a 

in  blood  at  lung  exit 

B 

in  brain 

CSF 

in  cerebrospinal  fluid 

I 

in  inspired  air 

P« 

PH 

a 

of  blood  at  lung  exit 

vB 

of  blood  at  brain  exit 

vT 

of  blood  at  tissue  ex'* 

CSF 

of  cerebrospinal  fluid 

Q 

blood  flow 

liters/min 

cardiac  output 

B 

cerebral 

N 

normal  (resting) 

AQ 

change  in  blood  flow 

liters/min 

— 

cardiac  output 

B 

cerebral 

(co2) 

due  to  carbon  dioxide 

(o2) 

due  to  oxygen 

r 

time  constant 

minutes 

1 

tor  cardiac  output 
response 

2 

for  cerebral  blood 
flow  response 

f 

time 

minutes 

Vi_ 

* 
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i  Subscripts 

Char¬ 

acter 

Loca¬ 

tion 

(i) 

Species 

0) 

Definitiotx 

Units 

T 

blood  transport  delay 

minutes 

aB 

from  lung  to  brain 

1 

aT 

from  lung  to  tissue 

I 

vB 

from  brain  to  lung 

vT 

from  tissue  to  lung 

ao 

from  lung  to  carotid 
body 

(1) 

in  appropriate  segment 
of  total  path 

V 

i 

gas  flow  rate  (alveolar) 

liters  (BTPS)/min 

E 

expiratory 

I 

inspiratory  (i.e., 
ventilation) 
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Appendix  B 

"NORMAL"  PARAMETER  VALUES* 


aC02 

m 

0.510** 

k 

Si 

0.00132 

a°2 

m 

0.0240 

K’ 

m 

795 

m 

0.0130 

*B 

M 

1.000 

B 

m 

760 

TP 

‘•CSF 

m 

0.100 

(BHC03)b 

m 

0.5470 

Kj. 

m 

39.00 

(bhco3)b 

m 

0.5850*** 

m 

3.00 

(bhco3)t 

m 

0.5850*** 

*®B(C02) 

m 

0.050 

(bhco3)CSf 

m 

0.5850 

**b(o2) 

m 

0.050 

Dco„ 

4m 

m 

81.99  x  10~7 

^(cop 

m 

0.182 

Do 

U2 

- 

4.361  x  10~7 

^(op 

m 

0.215 

dn2 

m 

2.524  x  10“7 

K  1 

6.000 

fi(co2) 

m 

0.0GC0 

qbn 

m 

0.750 

fi(02) 

m 

0.2100 

rl 

m 

0.100 

fi(n2) 

- 

0.7900 

r2 

m 

0,100 

(Hb) 

m 

0.2000 

Some  parameter  values  were  assigned  in  writing  the  system 
equations,  e.g.,  buffer  slope  and  Haldane  effect  [Eq.  (3.1;], 
vascular  segment  volumes  [Eqs.  (8. 1C)— (£.14) ],  controller  gains 
[Eqs.  (9. 1)~(9.4) ] ,  etc. 

## 

Although  equations  distinguish  between  a  values  for  bloc-’, 
brain,  tissue,  and  CSF,  they  were  all  set  equal  in  this  initial 
exploration. 

•ft--*  * 

These  values  now  appear  to  be  too  high  [37],  and  will 
require  revision  in  future  expirations. 
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